function DCM = data_fitting(analysis,vd,s)
% Meta-modelling of Bayes-optimal responses (Newton's method)
% FORMAT DCM = spm_meta_model(DCM)
%
% store estimates in DCM
%--------------------------------------------------------------------------
% DCM.M   - meta-model specification
%      M: [1 x m struct]   - hierarchical inference model (cf DEM.M)
%      G: [1 x s struct]   - generative process (for spm_ADEM or spm_ALAP)
%      U: [n x N double]   - n prior beliefs over N samples
%     pE: [1 x 1 struct]   - prior expectation of meta-model parameters
%     pC: [1 x 1 struct]   - prior variance of meta-model parameters
%
% DCM.xY  - data structure
%      y: [N x p double]   - N samples of a p-variate response
%     X0: [N x q double]   - q-variate confounds
%     dt: [1 x 1 double]   - size of time bin for each sample
%      Q: {[N x N double]} - precision component[s]
%
% DCM.xU  - input structure
%      u: [r x N double]   - r-variate input (hidden causes G in DEM)
%
% Computes (and stores in DCM_MM_???)
%--------------------------------------------------------------------------
% DCM.DEM - Inference (with MAP parameters)
% DCM.Ep  - conditional expectation
% DCM.Cp  - conditional covariances
% DCM.Eh  - conditional log-precision
% DCM.Ey  - conditional response
% DCM.F   - log-evidence
%
% This routine illustrates Bayesian meta modelling - the Bayesian inversion
% of a model of a Bayesian observer. This requires the specification of two
% models: An inference model used by the subject (specified by a DEM
% structure) and a meta-model (specified by a DCM structure). The inference
% model is completed by a response model to furnish the meta-model; where
% the response model takes the output of the (active) inference scheme
% specified by the DEM and generates an observed (behavioural or
% neurophysiological) response. Crucially either the inference model or
% the response model or both can have free parameters - that are optimised
% using Bayesian nonlinear system identification in the usual way.
%
% Although this routine is a function, it is expected that people will fill
% in the model-specific parts in a local copy, before running it.  The
% current example uses a model of slow pursuit and generates synthetic data
% (responses) to illustrate how it works. To replace these simulated data
% with real data, simply specify the DCM.xY (and xU fields) with
% empirical values. If other fields do not exist, exemplar fields will be filled in.
%
% The conditional density of the parameters and F values (log-evidence) can
% be used in the usual way for inference on parameters or Bayesian model
% comparison (as for other DCMs)
%__________________________________________________________________________
% Copyright Eduardo Aponte (C) 2008 Wellcome Trust Centre for Neuroimaging



spm_figure('GetWin','Experimental Data');
clf
subplot(2,1,1);
analysis.plotTrial();

%%
if nargin == 1
    [vd,s,ox,ow] = analysis.normalizedData();
elseif nargin > 1
    [ans,ans,ox,ow] = analysis.normalizedData();
end

% Transformed to degrees

v2d = @(x) sign(x).*acosd(16./(4 * sqrt(16 + x.^2)));

vd = v2d(vd);
s  = v2d(s);
ox = v2d(ox);
ow = v2d(ow);

% Normalize to unit
% =========================================================================

nc = 15;

vd = vd/nc;
s = s/nc;
ox = ox/nc;
ow = ow/nc;
c = -8/nc:1/nc:8/nc;

% INFERENCE MODEL: DEM.M and DEM.G
%==========================================================================

% create illustrative (M,G) if not specified (EXAMPLE)
%----------------------------------------------------------------------

% hidden causes and states
%======================================================================
% x    - hidden states:
%   x.o(1) - oculomotor angle
%   x.o(2) - oculomotor velocity
%   x.x(1) - target angle - extrinsic coordinates
%
% v    - causal states: force on target
%
% g    - sensations:
%   g(1) - oculomotor angle (proprioception)
%   g(2) - oculomotor velocity
%   g(:) - visual input - intrinsic coordinates
%----------------------------------------------------------------------

%%
% Set-up
%======================================================================
M(1).E.s = 1/2;                               % smoothness
M(1).E.n = 4;                                 % order of
M(1).E.d = 1;                                 % generalised motion


% sensory mappings with and without occlusion
%----------------------------------------------------------------------


% Generative model (M) (sinusoidal movement)
%======================================================================
% Endow the model with internal dynamics (a simple oscillator) so that is
% recognises and remembers the trajectory to anticipate jumps in rectified
% sinusoidal motion.

% slow pursuit following with (second order) generative model
%----------------------------------------------------------------------
x.o = [vd(1);0];                                  % motor angle & velocity
x.x = vd(1);                                      % target location

% level 1: Displacement dynamics and mapping to sensory/proprioception
%----------------------------------------------------------------------
elasticity = 10;
tmpf = @(x,v) exp(-(c' - (x.x - x.o(1))).^2);
tmpdf = @(x) -2*(c' - (x.x - x.o(1)));

M(1).f = @(x,v,p) [x.o(2); (v - x.o(1))/4 - x.o(2)/2; v - x.x];

M(1).g  = @(x,v,p) [x.o(1);x.o(2); tmpf(x).*~(x.x > (ox-abs(ow)/2) && x.x < (ox+abs(ow)/2))];

M(1).gx = @(x,v,p) [1,0,tmpf(x)'.*tmpdf(x)'.*~(x.x > (ox-abs(ow)/2) && x.x < (ox+abs(ow)/2))
    0,1,zeros(1,17)
    0,0,-tmpf(x)'.*tmpdf(x)'.*~(x.x > (ox-abs(ow)/2) && x.x < (ox+abs(ow)/2))]';
M(1).gv = @(x,v,p) sparse(19,1);

M(1).fx = @(x,v,p) [0,-0.25,0; 1,-0.5,0;0,0,-1]';
M(1).fv = @(x,v,p) [0;0.25;1];


M(1).x = x;                                   % hidden states
M(1).V = exp(4);                              % error precision
M(1).W = exp(4);                              % error precision


% level 2: With hidden (memory) states
%----------------------------------------------------------------------
M(2).f  = @(x,v,p) [x(2); -x(1)]*v;
M(2).g  = @(x,v,p) x(1);

M(2).fx = @(x,v,p) [0,1;-1,0]*v;
M(2).fv = @(x,v,p) [x(2); -x(1)];

M(2).gx = @(x,v,P) sparse([1,0]);
M(2).gv = @(x,v,P) sparse(1,1);

M(2).x  = [0; 0];                             % hidden states
M(2).V  = exp(4);                             % error precision
M(2).W  = exp(4);                             % error precision

% level 3: Encoding frequency of memory states (U)
%----------------------------------------------------------------------
M(3).v = 0;
M(3).V = exp(4);


% generative model (G)
%======================================================================

% first level
%----------------------------------------------------------------------
damp = 10;

%elasticity = 0;
G(1).f = @(x,v,a,p) [x.o(2); a - x.o(2)/damp - x.o(1)/elasticity; v - x.x];
G(1).g = @(x,v,a,p) M(1).g(x,v,p);

G(1).fx = @(x,v,a,p) [0,-1/elasticity,0;1,-1/damp,0;0,0,-1]';
G(1).fv = @(x,v,a,p) [0, 0, 1]';
G(1).fa = @(x,v,a,p) [0, 1, 0]';

G(1).gx = @(x,v,a,p) M(1).gx(x,v,p);
G(1).ga = @(x,v,a,p) sparse(19,1);
G(1).gv = @(x,v,a,p) sparse(19,1);

G(1).x = x;                                  % hidden states
G(1).U = sparse(1,[1 2],[1 1],1,19)*exp(0);  % motor gain

% second level
%----------------------------------------------------------------------
G(2).v = 0;                                  % exogenous force
G(2).a = 0;                                  % action force

%======================================================================

% Check generative model
%--------------------------------------------------------------------------
M      = spm_DEM_M_set(M);


% Experimental input (C)
%==========================================================================

xU.dt = 16.6666; % ms
N     = numel(vd);
xU.u  = vd;


% Priors (U)
%==========================================================================


% (EXAMPLE) prior beliefs about target frequency (w)
%----------------------------------------------------------------------
U = zeros(1,N) + 2*pi/analysis.FpC(1);

% Data and confounds
%==========================================================================

% Generate simulated data (EXAMPLE)
%----------------------------------------------------------------------

% meta-model and true parameters
%----------------------------------------------------------------------
MM.M    = M;
MM.G    = G;
MM.U    = U;


xY.y    = s';


% DCT confounds
%==========================================================================
%xY.X0  = spm_dctmtx(N,1);

% and [serial] correlations (precision components) AR model
%--------------------------------------------------------------------------
xY.Q   = {spm_Q(1/2,N,1)};


% META-MODEL parameters
%==========================================================================
% These parameterise the inference scheme function (IS) at the bottom of
% this script - this functions defines how the parameters are used and
% therefore optimised. Change this function to specify the meta-model


% prior expectations (EXAMPLE)
%----------------------------------------------------------------------
%pE.W  = -0.6708;
%pE.W  = ;
pE.W  = 2;
%pE.a  = 0.6137;


% prior covariance (EXAMPLE)
%----------------------------------------------------------------------
pC.W  = 1;
%pC.a  = 1;


% This completes the meta-model specification. The fields are now assembled
% and passed to spm_nlsi_GN (nonlinear system identification using Gauss-
% Newton-like gradient ascent).
%==========================================================================


% META-MODEL (MM)
%==========================================================================
MM.M   = M;
MM.G   = G;

% hyperpriors (assuming about 99% signal to noise)
%--------------------------------------------------------------------------
hE     = 8 - log(var(spm_vec(xY.y)));
hC     = exp(-4);

% Meta-model
%--------------------------------------------------------------------------
MM.IS  = @IS;
MM.pE  = pE;
MM.pC  = pC;
MM.hE  = hE;
MM.hC  = hC;

%%
% Model inversion
%==========================================================================

[Ey,DEM] = IS(pE,MM,xU);
spm_figure('GetWin','Generated Data');
spm_DEM_qU(DEM.qU)
keyboard

[Ep,Cp,Eh,F] = spm_nlsi_GN(MM,xU,xY);

[Ey,DEM] = IS(Ep,MM,xU);

% integrate (A)DEM scheme with MAP etimates
%--------------------------------------------------------------------------

spm_figure('GetWin','Experimental Data');
subplot(2,1,2);
hold on
hs = plot((0:numel(s)-1)/analysis.scn.refresh_rate,Ey,'r');
hv = plot((0:numel(vd)-1)/analysis.scn.refresh_rate,s','k');
legend([hs,hv],'Eye trace','Stimuli','location','NorthEast');
title('Data: Unit amplitude')
ylabel('Horizontal axis');
xlabel('Secs.');

% store estimates in DCM
%--------------------------------------------------------------------------
DCM.DEM = DEM;                  % Inference scheme (with MAP parameters)
DCM.M   = MM;                   % meta-model
DCM.xY  = xY;                   % data structure
DCM.xU  = xU;                   % input structure
DCM.Ep  = Ep;                   % conditional expectation
DCM.Cp  = Cp;                   % conditional covariances
DCM.Eh  = Eh;                   % conditional log-precision
DCM.Ey  = Ey;                   % conditional response
DCM.F   = F;                    % log-evidence

return


% inference scheme x = IS(p,M,U)
%==========================================================================
function [y,DEM] = IS(P,M,U)

% parameterise inference model (EXAMPLE)
%--------------------------------------------------------------------------

M.M(2).W = exp(P.W);
%M.G(1).U = sparse(1,[1 2],[1 1],1,19)*exp(P.a);  % motor gain

% reset random number generator (for action schemes)
%--------------------------------------------------------------------------
% Bayesian inference (spm_DEM, spm_LAP, spm_ADEM or spm_ALAP)
%--------------------------------------------------------------------------
DEM.M  = M.M;
DEM.G  = M.G;
DEM.U  = M.U;
DEM.C  = U.u;
DEM    = spm_ADEM(DEM);

% (Hidden) behavioural or electrophysiological response (EXAMPLE)
%--------------------------------------------------------------------------
x      = DEM.Y(1,:);
y      = x';

return